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A multicomponent extension of our recent theory of simple fluids [ U.M.B. Marconi and S. Mel- 
chionna, Journal of Chemical Physics, 131, 014105 (2009) ] is proposed to describe miscible and 
immiscible liquid mixtures under inhomogeneous, non steady conditions typical of confined fluid 
flows. We first derive from a microscopic level the evolution equations of the phase space distribu- 
tion function of each component in terms of a set of self consistent fields, representing both body 
forces and viscous forces (forces dependent on the density distributions in the fluid and on the 
velocity distributions). Secondly, we solve numerically the resulting governing equations by means 
of the Lattice Boltzmann method whose implementation contains novel features with respect to 
existing approaches. Our model incorporates hydrodynamic flow, diffusion, surface tension, and the 
possibility for global and local viscosity variations. We validate our model by studying the bulk 
viscosity dependence of the mixture on concentration, packing fraction and size ratio. Finally we 
consider inhomogeneous systems and study the dynamics of mixtures in slits of molecular thickness 
and relate structural and flow properties. 

PACS numbers: 47.11.-j, 47.61.-k, 61.20.-p 



I. INTRODUCTION 



In recent years, the experimental study of the flow of liquids in nanometric structures has become an important 
and rapidly evolving discipline driven by technological applications, such as lab-on-a-chip systems, electrophoresis and 
electro-osmotic pumping. The mixing and separation of different molecules according to their size or chemical and 
transport properties is a major issue in chromatography, nanofluidic logic gates and electrokinetics. Many applications 
are foreseen also in biotechnologies, because very small quantities of liquid components arc sufficient for analysis and 
synthesis, including the possibility of isolating and analyzing or manipulating a single molecule in a single nanochannel 
or surface with nanoscale features [l|, . 

These technological advances require a better theoretical understanding of the fluid properties in confined geometries 
and under non-equilibrium conditions. At the nanoscale, new areas of physics, chemistry and materials science 
come in. The systems become more surface-like and molecules never explore a bulk-like environment and because 
inhomogeneities have a great impact some interesting phenomena arise. Whereas a large body of information has been 
accumulated in the last thirty years concerning the physics of inhomogeneous fluids under thermodynamic equilibrium 
conditions [||, the state of the art of flowing fluids is not so advanced 0-0] ■ 

It is crucial to understand how the interplay between the various forces occurs and how it changes when going 
from macrosystcms to nanosystems. It is clear from dimensional considerations that the relative importance of forces 
changes with the typical dimensions of the systems and we expect that the smaller the sample under scrutiny the 
more important the role of surface forces is. Not only nanoscale flows are dominated by viscous effects if turbulence is 
absent, but wherever frequent molecule-channel wall collisions occur besides molecule-molecule collisions, we expect 
a modified frictional resistance Q. 

The generalization of fluid dynamics from pure to multicomponent fluid mixtures composed of different components 
or species requires the introduction of some new concepts. While traditional computational fluid methods involve 
some ad hoc extrapolation of Navier-Stokcs equation to confining geometries, our approach is based on a microscopic 
kinetic equation incorporating the effects of the inhomogeneities in its structure. 

The major difference with respect to one-component fluids is the phenomenon of concentration diffusion. We shall 
study a multicomponent kinetic equation |9l4l3j that has been considered in the past by several authors with the 
specific goal of computing the transport coefficients from a microscopic approach. Many of these studies, which were 
based on the Boltzmann equation or on its extension to the dense case, the Boltzmann-Enskog equation, require some 
analytical and numerical effort, so that a simpler treatment of the interactions has been proposed. A very popular 
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approximation is represented by the phenomenological Bhatnagar-Gross-Krook (BGK) equation whose merit is to 
reduce the complexity of the original Boltzmann collision kernel f2 Q ^ between species a and (3 by introducing a simple 
relaxation-time ansatz The BGK, in spite of being even less accurate than the original Boltzmann equation, 

has enjoyed a great popularity especially in conjunction with the Lattice Boltzmann method (LBM) thanks to the 
simplicity of the collision kernel leading to a considerable speed-up of the numerics [l9j . 

However, in the multicomponent case the choice of the form of the BGK relaxation term is not unique and hard to 
infer from the original f2 Q/3 , and the literature abounds of several proposals poM23 | . 

In order to describe with sufficient accuracy the fluid structure at length scales comparable with the size of the 
particles we shall resort to methods similar to those of density functional theory (DFT) employed in the study of 
equilibrium and non equilibrium properties [25l - |28j . In the case of hard-core fluids, DFT and its dynamical extension 
gives excellent results and can be extended to more realistic fluids by using the van der Waals picture of decomposing 
the total inter-particle potential into a short-range repulsive potential and a long-range attractive potential tail. The 
first is treated by means of a reference hard-sphere system whilst the second is considered within the random phase 
approximation (RPA). Such a decomposition can describe phase separating systems, wetting phenomena, two phase 
interfaces, but predicts transport coefficients of purely hard-core systems [29(. 

The hard-core part of the interaction is treated within the revised Enskog theory (RET) [3(j that considers the 
non-local character of the momentum and energy transfer and takes into account the static spatial correlations among 
particles. This level of theory does not incorporate the time correlations which are responsible for memory effects and 
hydrodynamic contributions to the transport coefficients {3l| . 

With respect to the dynamical DFT (32, HH , the present theory preserves Galilei invariance and therefore displays 
full hydrodynamic behavior in the limit of slowly varying fluctuations about the equilibrium reference state. 

With the aim of deriving a practical numerical approximation to study mixtures in inhomogeneous situations, we 
extended the method of Dufty and coworkers |34j to the multicomponent case. The method is a compromise between 
the RET, of which it retains the accuracy as far as the momentum and energy transfer are involved, and the much 
simpler BGK, which we employ to evolve the non-hydrodynamic moments of the distribution functions (3^413 ■ 

The final product of our theory is a coupled system of simplified equations for the density distributions of individual 
species describing both the streaming and the collisional stages. These equations are solved by an appropriate extension 
of the LBM algorithm to include both particle-particle interactions and particle- wall interactions. A simple analysis of 
the equations is used to derive explicit expressions both for equilibrium thermodynamic quantities, such as pressure, 
compressibility, etc., and for non equilibrium transport coefficients. It is important to stress the difference between the 
present work and other LBM based approaches for non-ideal fluids, where only attractive interactions are accounted 
for by the so-called Shan-Chen pseudo-potential, whose justification is purely mcsoscopic, while transport coefficients 
enter as free parameters of the theory [H, |39[ . 

Interestingly, over two decades ago Davis and coworkers proposed on phenomenological grounds a local average 
density model (LADM) of viscosity and diffusivity, based on the principle of local averaging of the density over the 
volume of a hard sphere |4(j. The present approach derives similar expressions for the transport coefficients, but differs 
from the LADM in the way it solves the resulting evolution equations. Pozhar gave a more rigorous foundation to 
the transport theory in non-uniform systems, but the resulting theory could not be developed into a computationally 
tractable method [4l|, |42[. In the present work, we claim that our microscopically founded theory lends itself to 
relatively simple numerical solutions based on the use of the LBM even in the presence of complex geometries. The 
numerical application regards fluids confined to a spacing of a few monolayers between solid surfaces, where they may 
behave very differently from bulk conditions and exhibit unusual properties such as a dramatic increase in viscosity. 

The paper is organized as follows: in Sec. |TI]we introduce the model and the evolution equations for the individual 
distribution functions of the multicomponent system, and in Sec. IIIII we discuss the approximations involved in the 
theory. In Sec. IIVI we test the predictions of the resulting evolution equations using a simple confining geometry. 
In Sec. [V] we validate the resulting equations numerically and simulate the flow of a mixture in narrow channels. 
Finally in Sec. [VI] we present our conclusions and perspectives. Two appendices are added to the paper in order to 
illustrate the relationship between the collision kernel and the non-ideal part of the chemical potential and discuss 
some technical details related to the discretization procedure. 



II. KINETIC EQUATIONS FOR THE DISTRIBUTION FUNCTIONS 

We consider a fluid mixture of N a particles of mass m a , with a — 1,M, where M represents the number of 
components, mutually interacting with the potentials U a "(r, r'). The probability density of finding a particle of type 
a at position r with velocity v at time t is proportional to the the single-particle distribution functions, / a (r,v,t), 
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normalized as follows: 

N a = J dv J dvf a (r,v,t). (1) 

The evolution of f a is given, in principle, by the exact BBGKY hierarchy of dynamical equations for the reduced 
distribution functions [9j, whose first member is the following kinetic equation for the a component: 

d t f a (r, v, t) + v ■ V/ Q (r, v, f) + ^ ■ (r, v, t) = £ ^(r, v, i), 

(2) 

where F a (r) is an external velocity independent force field acting on component a and fl al3 represents the effect of 
the interactions among the fluid particles of type a and /?. In principle, the exact form of fi"' 3 is known, but it 
involves the two-particle distribution functions, which in turn depend on the three particle distribution functions. In 
order to close this hierarchy different approximations have been proposed, the first historically being the Boltzmann 
equation expressing il a ^ in terms of products of single-particle distributions. Other specific choices of f2 Q ' 9 in eq.@ 
correspond to different kinds of approximations and lead to equations such as the BGK equation [l4| , the Vlasov or 
RPA equation Q and the RET equation [3p( ] . 



A. Macroscopic variables. 

We seek a description in terms of a small number of macroscopic variables which arc functions of the fixed position 
r and the time t. The procedure leading from eq.© to the evolution equations in terms of these variables, closely 
follows the analysis of the one-component case |36[ and allows to derive a set of hydrodynamic equations, in the case 
of a M-componcnt fluid mixture. In D spatial dimensions, there exist L = [D + M + 1) conservation laws, for the 
mass of each species, global momentum and energy. Thus, one can consider the evolution of L linear combinations of 
velocity moments of f a which arc identified with the hydrodynamic fields of the mixture. 

In the following, we shall adopt the Einstein summation convention for repeated Cartesian indexes, while sums over 
components are written explicitly. Let us first write the statistical expressions for the mass density of component a: 



p a {v, t) = m a n a (r, t) = m a J dv/ Q (r, v, i), (3) 

where n a is the associated number density. The total average mass density is given by: 

p(r,t)=Y / P a (r,t), (4) 

a 

while the total number density is n = ^ a n a . The first velocity moment of the phase-space distribution corresponds 
to the average local velocity, u Q , of the a component 

uQ(r '' ) = ^)/ dVV/Q(r ' V '* ) - (5) 
We also consider the local barycentric velocity: 

u r,t = V ' i— 1— i 6 

and the so called dissipativc diffusion current, measuring the drift component a with respect to the center of mass 
velocity: 

J Q M) = p a (r,t)[(u a (r,t) - u(r,t)] = p a (r, i)w Q (r, t) (7) 
and having the property J2 a J a ( r ; *) = 0. 
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B. Balance equations. 

The macroscopic variables p a ,u a ,u obey a set of evolution equations derived from eq. ([2]) as follows. Upon 
multiplying by m a and integrating over the velocity degrees of freedom both sides of eq. ([2} one obtains the continuity 
equation for the mass density of the species a: 

d t p a {v, t) = -V • (>(r, t)u a {v, t)) = -V • (>(r, t)u(r, t) + J Q (r, tj) . (8) 

After summing over components, one obtains the global mass continuity equation: 

d t p(v,t) + W- (p(r,t)u(r,f))=0. (9) 

In establishing eq. ([HJ we have used the property that the interaction does not affect the partial number densities, i.e. 
/ dvW 13 = 0. The momenta of the single species are not locally conserved quantities, due to the interaction between 
different components. To derive the associated governing equations we first multiply cq. ^) by m Q (v — u) and then 
integrate over the velocity with the result: 

dt [p a (r , i)< (r, *)] + Vi (p a (r , t)uf (r , t)uf (r, t) - p Q (r , t)< (r , tW? (r , t) 



+ ^£W*) +E c f ( r >*)- ( 10 ) 

In analogy with pure fluids, in eq. (|10[) we have introduced the kinetic contribution to the partial stress tensor: 

7rg(r,i) = p Q (r,<) < («i - Ui)(vj - u^-) >= ra a / dw{ Vl - u i )(« j - i^)/ Q (r, v, i). (11) 



Notice that in order to ensure that 7r£* is independent of a superimposed uniform translation of the entire flowing 
fluid its definition contains the difference {vf — itj). In eq. pU|) we have also defined the projection of the collision 
term onto the velocity: 

C a,3 (r, t)=m a J dvvVL al3 {r, v, t), (12) 

representing the contribution to the stress tensor arising from interactions. This integral will be explicitly computed 
in Sec. IIII Al for a specific model. Finally, the term p a w?w° l in the l.h.s. of eq. (fT0|). has nothing to do with 
viscous effects, but represents an additional stress of dynamical origin due to the different average velocities of the 
two components |44| and vanishes in the pure fluid. Summing over species eq. (|10[) we obtain the local conservation 
law for the total momentum, pu, 

dt[p{r,t) U] (r,t)] + Vi (p(r,t)ui(r,t)uj(r,t)]) 

a a/3 

with the total kinetic component of the stress tensor given by: 

^(MHE^-M)- ( 14 ) 

a 

The set of eqs. © and (fT5]) is not determined unless additional relations specifying the pressure tensor (fT4")) in terms of 
the selected macroscopic fields arc given. Hydrodynamics at the Navier-Stokes level closes the equations by assuming 
a simple linear relation between fluxes and derivatives of the fields via phenomenological transport coefficients (such 
as the viscosity and conductivity) whose values are not known a priori, but which must be taken from experiments. 
Kinetic theory, working at the level of the non-equilibrium distribution functions, avoids such assumptions. 



III. SIMPLIFIED KINETIC MODEL 



In a nutshell, our strategy consists in solving the kinetic equations rather than the balance equations ((5J) and (fT5]) . It 
involves concepts borrowed from the lattice Boltzmann method (LBM) and from density functional theory (DFT). The 
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LBM can be regarded as an efficient approach to treat inhomogencous systems and multiphase flows and provides 
an alternative to the solution of the Navier-Stokes equations. Conditions at solid boundaries are relatively easily 
implemented using the LBM, making it well suited for handling complex geometries. The LBM becomes particularly 
useful to treat non-local effects, such as wetting layers, interfaces, layering etc., if used in conjunction with a simplified 
kinetic approach suggested by Dufty and coworkers in the case of one-component fluids [34j . 

To introduce the simplified kinetic model for mixtures, we refer to ref. (36j . where we discussed the simpler one- 
component case and reduced the full transport equation ([2J to a simpler form. The extension to the multicomponcnt 
case requires some subtleties which are needed in order to satisfy a key requirement of this approximation: the 
modified kinetic equation and the original one must share the same equations for the partial densities and momenta, 
i.e. the distributions must have the same low order projections. Doing so, we are able to treat exactly the contribution 
of the collision operator to the hydrodynamic equations. As a result, we need to make approximations only on that 
part of the collision operator which acts on the higher order, non-conserved modes. This latter contribution can be 
handled in the spirit of some cruder approximation, such as the BGK-like approximation [l8j . Since the velocity 
distributions in the final equilibrium state are known it is easy to construct the BGK collision operator in such a 
way that it verifies some basic properties of the original problem, namely particle and momentum conservation. Our 
two-component extension of the simplified kinetic equation reads: 

d t f a (r, v, t) + v ■ V/ Q (r, v, t) + • A/« ( r , v , i) 



1 V"(r,v,i) 



(v - u) • C ap (r, t) - oj[f a (v , v, t) - (r, v, t)\ , 



k B T n a (r,t) 

(15) 

where w is a collision frequency which depends on the two species and 



and 

ih a fr v A = ih a frvtlllJ i K —L-L 



(17) 

The proposed form (|15p differs from the existing literature about mixtures in a few aspects. In the first place, the 
mixtures were treated cither by including the full collision operator or by assuming the BGK to be a valid replacement. 
The simplified kinetic method, to the best of our knowledge, was never applied to mixtures. Our choice for the BGK 
term is based on a single relaxation time for the two components and represents perhaps the minimal model satisfying 
the conservation laws compatible with the Dufty and coworkers strategy. 

A few comments are in order. First of all, eqs. ©, (JTUJ) and (|T5)) arc immediately recovered after multiplying cq. (fT5)) 
by m",m°(v-u) and integrating over the velocity, v. Second, the coupling between the evolution equations of the two 
distributions f a occurs through the nonlinear collision kernel and this fact is also reflected in the BGK representation 
of that kernel. Equation ([2]) contains two terms 51 Q ' 3 , one for each type of interaction, but in the modeling of the BGK 
relaxation time approximation it is convenient to choose a single collision term. This assumption is valid not only 
because the coupling between the two distributions is achieved via the barycentric velocity u contained in %j) a , but 
also through the term C Q ^. Third, the classical BGK approximation when applied to a one-component fluid would 
involve the difference between the distribution f a and the local Maxwcllian ip a . In the present case a non exponential 
prcfactor in eq. (fTTj) has been introduced and serves as to "orthogonalizc" the term tu[f a — ijjj] to the "true" collisional 
terms, C™^, so that the BGK does not lead to modifications of the hydrodynamic equations. Also notice that the first 
velocity moment of ip a is n Q u, whereas the first moment corresponding to is n a u a . The present ansatz reduces 
to our one-component equation [36j | in the limit i/'" ~~ ^ V 1 "- The above choice satisfies the indifferentiability principle 
which dictates that when all physical properties of the species are identical, the total distribution / = f A + f B obeys 
the single species transport equation. 

A. Interactions 

We consider the interaction term featuring in eq. ([2]) which according to the BBGKY hierarchy can be written as 
fi°*(r, v, t) = -Lv v -Jdv> J dr'VrU^iv - v')ff{r, r', v, v', i), (18) 
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where /^(r, r', v, v', t) represents the two-particle distribution function and the U a ^((r — r') are the pair potentials 
between particles of species a and f). 

In general, the intcrmolecular potential, U al3 , contains both a repulsive and an attractive contribution. The latter 
plays a major role in determining the thermodynamic properties of the system. In the following we shall separate 
the repulsive and the attractive contributions by representing the short range repulsion between the particles through 
repulsive hard sphere potentials, characterized by the three diameters ctaa, ctbb and gab = {&AA + obb)/2 plus an 
attractive contribution: 

fi<*(r,v,t) = 0^(r,v,t) +fCi(r,v,t). (19) 
The repulsive part is treated within the following RET approximation |45j : 

n?£ p (r,v a y,t) =oipf dv? J dk6(k-v a/3 )(k-v^)x 
{g a p(r, r - a a0 k, t)f a (r, v a ,t)f?(r - a^k, v , t) 

-g a p(r, r + a a pk, t)f a (r, v a ,t)f(r + a a/3 k, , t)}, (20) 
where v Q( g = (v Q — v* 3 ), while v a and v' 3 are scattered velocities given by 

2mP 



mr + m 



(k- v Q/3 )k 



^=v + -^— ,(k-v Q/3 )k. ( 21 ) 



The quantities g a p(r,r ± er a| gk) are the three hard sphere pair correlation functions evaluated when the particles of 
species a and j3 are at contact. The RET collisional term describes the large momentum binary exchanges due to close 
encounters between the cores, whereas the attractive term can be seen as a slowly varying average force, involving 
small momentum exchanges. In detail, the attractive tails of the potentials are treated in a mean field approximation, 
rewriting the second term of the r.h.s. of eq. (|19j) by factorizing the two-particle distribution functions 

ff(r, v, r', V, t) » <Mr, r', t)f (r, v, t)f(r', v', t), (22) 
that is, by neglecting velocity correlations completely: 

n^r.v.t) « /rfv'| rfr'V r t/ a ^((r-r') 5a Mr,r',i)r(r,v,O/ (r',v',i), (23) 

where g a p{r, r', £) is the inhomogeneous the pair correlation function of the reference hard-sphere system. Integrating 
w.r.t. the velocity, we obtain the RPA approximation for the attractive term: 

0^ r (r,v,i) = _^!!M . V t ,r(r,v,i), (24) 

where we introduced the molecular fields 

G<*(r,t) = -J d/^(r',t) ffQ/3 (r,r')V r t/ a Q t i(r-r'). (25) 

In order to evaluate the interaction contribution to the momentum balance equation, we use the definition (|12[) 
together with eqs. (|20|) and ([24]) . After some algebra, using explicitly the collision rules eqs. (|2"Tj) we obtain the 
following expression: 

Cf{v, t) = -Zucpffip J rfv Q J dv* 9 J dkQ(k • v a/3 )(k • v Q/3 )(k • v a/J )k, 

g a p(r, r + a Q(9 k, t)f a (r, v a , t)/^(r + a a/3 k, v* 9 , i) + Gf (r, t)n a (r, t) (26) 
where we have introduced the reduced mass fi a p = ^J^p ■ To proceed further, we invoke a local equilibrium 



approximation and replace the distribution functions within the integrals (|26j) by local Maxwellians 

r (r ,v, f) K r (r ,v, t) s „ nM)[ ^|^exp(- !rgi__gM))! ) (27) 
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and perform the integrals involved in eq. ([20]) . by expanding u a (r,t) about u(r,t). We also expand the functions at 
point (r + dap) about their values at r up to first order in the differences. The procedure is completely analogous to 
the one followed in ref. (3(| and leads to the final result : 

C Q,3 (r, t) = -a 2 aj3 j dkkg al3 (r,r + a a pk,t)n a (r,t)n f3 (r + a af3 k,t) x 



[k B T _ 2 J 2 ^ BT k ■ [u B (r + a a0 k, t) - u«(r, t)]) +G^(r, t)n a (r, t). 



(28) 



B. Decomposition of the interaction term into specific force contributions. 

It is instructive to rewrite the term featuring in formula (|28[) as a sum of more specific forces: 

J2 C a/3 (r, t) = n Q (r, t) (v a ' m f (r, t) + F Q ' dra9 (r, t) + F a > mscous {r, tj) . (29) 

8 

We identify the force, F a,m f , acting on the a particle at r due to the influence of all remaining particles, being the 
gradient of the so-called potential of mean force: 

F<*' m f(T,t) = -k B Tj2<r 2 «0 I dkkg a p(r,T + a aB k,t)n p {v + a aB k,t) + ^G a/J (r,t) (30) 

d 9 

a drag force with local character: 

Ft draa (r,t) = -£ 7 |f (r)[<(r) - u J( r )] (31) 

8 

and a viscous force of non local character: 

K MSCOU S (r,i) = £ / dr'^(r,r')[<(r')-<(r)]. (32) 
8 J 

We have further defined a friction tensor: 

7 ^(r) = J dv'H* 8 (r,r>) (33) 

and a viscosity kernel: 



H tj P {r,r) = 2 ] J g a p{r, r , t)n p (r , t)5{\r - T \-<r aft ) ' ( 34 ) 

whose relation with the transport coefficients will become clear in the following. 

The term F°'™^ is related to the gradient of the excess chemical potential of species a, nf nt (r), due to the interactions 
among the particles. In Appendix A we show that in the case of slowly varying density profiles one has the result 

F^(r,i) = -V/Ct(r,*). (35) 

The drag term is proportional to the velocity difference of the two species and is present even in the absence of 
velocity gradients. The viscous term, on the other hand, is proportional to velocity gradients. We can consider the 
case u A = const and u B = const and u A ^ u B and uniform densities so that the tensor 7?. becomes isotropic 

F A ^™ 9 (r, t) = - 7 AS (u A - u s ). (36) 

By using the definition (|33)) we obtain 

j AB = % -{^p, AB k B TYl 2 o\ B g AB n B . (37) 
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It is interesting to show how the microscopic viscosity kernel (|34)) is related to the viscosity coefficient. Let us 
assume a simple case u A = u B = u with 

u(r,t) = (0,0, «*(£,?/)) (38) 

and uniform densities, n A and n B . In this case the only non vanishing component of the viscous force is directed 
along z. After some lengthy, but straightforward algebra (see ref [HI) we obtain: 

K „ = 4. J2 V2x^fc fl T^ Sofl + ^) (39) 

13 y 

which allows to identify the collisional contribution to the shear viscosity coefficient by using the macroscopic relation 

=. (C) (^ + ^) (40) 

a 

from which we can write the explicit formula: 



4 

15 

a/3 ufl 



V (C) = E ^3 = ^ E y/^c.pk B Tn a nPai p g a p. (41) 



In the case of pure fluids, this expression of the viscosity reduces to the one predicted by Longuet-Higgins [46] , which 
is known to describe successfully fluids of hard molecules even when the radial distribution function is momentarily 
isotropic. In spite of the fact that the trial distributions correspond the the local equilibrium state, the singular nature 
of the interactions allows a finite flux of momentum. Notice that the RPA attractive term does not contribute to the 
transport coefficients, but only to the potential of mean force. 



C. Fischer-Methfessel prescription for the pair correlation functions of a non uniform system. 

We have deduced the formulas for the viscosity and drag coefficient in the case of systems of uniform densities. In 
order to apply our theory to non-uniform cases we need to specify how the pair correlation functions are modified 
with respect to the bulk. The pair distribution function g a p appearing in the above formulas is constructed according 
to a generalization of the prescription put forward by Fischer and Methfessel (FM) (4?J . One first defines the smeared 
densities n a {r) via a uniform averaging of n a (r) over spheres of volume V a = ir<J% a /& centered at the point R Q( g 
located between the centers of the two spheres at r Q and rp, respectively: 

IW = r -^- (42) 

The smeared densities are [48j : 

n«(R a p) = ^- [ drn a (r-R a p). (43) 

"a JV a 

The contact values of the pair correlation functions arc computed using the bulk expressions obtained from an extension 
[49l [50I ] of the Carnahan-Starling equation to mixtures: 



1-& 2 a a p (1 — 43) 2\ cr Q(3 J (l-£ 3 ) d 

where the functions £„ are evaluated for the values of the smeared densities at point R a ^ between the two particles 
at contact: 

e n (R^) = ^E^( R ^K- ( 45 ) 

a 

The pair correlations at contact are evaluated in the inhomogeneous system using the generalized FM ansatz: 

g a p(r a ,rp; \r a - r fj \ = a a p) = g b ^ k ({f„(Rap)». (46) 
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From the knowledge of g a p, one can immediately obtain the equation of state for the mixture: 



P 2 



nJ2n a nPg b ^ k al . (47) 



k B T 3 

After having made the connection between the collisional term C Q ^(r, t) and the three types of forces acting in the 
fluid eq. (p?9")) . we can use the large body of knowledge accumulated over the last twenty years about entropic forces in 
multicomponent systems (5l| . These forces, for instance, determine an effective attraction towards a flat wall of the 
particles with the larger radius induced by the presence of the small particles. For similar entropic reasons two large 
particles immersed in a sea of small particles experience an effective attraction. To the best of our knowledge, the 
interplay between entropic forces and viscous forces is not very well known and is worth to be explored. As an order 
of magnitude estimate entropic forced are typically ksTna AB , whereas the ratio between viscous forces and entropic 
forces is m u/v^ with = yJksT ' / ' p a p. 

Although recent versions of equilibrium DFT are more accurate than the FM method, it is more convenient to 
work within the FM approximation for a series of reasons: a) the FM approximation is very simple as compared to 



the Roscnfcld [52| prescription for pair correlations; b) our theory not only requires the effective potential V/if nt , but 
also the viscous and factional forces which are not given by the DFT method; c) computationally the FM method is 
much faster, although it becomes numerically unstable at high packing fractions. 



D. Numerical validation 



We consider a binary mixture of HS with <jaa — 4, <jbb = 8 and composition Xa = n A /(n A + n B ). The first 
validation of the numerical algorithm, succinctly described in Appendix B, is obtained by calculating the shear 
viscosity of the binary system and initially subjected to a sinusoidal shear wave of small amplitude in bulk conditions. 

The measured decay time is found to have the same temporal decay for both species, and the resulting viscosity 
curves are reported in Fig. [TJ The data are compared with the analytical expression for the decay time obtained 
by solving the linearized hydrodynamic equations and using the theoretical viscosity of eq. (|41|) . Notice that the 
viscosity increases monotonically with the packing fraction £ 3 of the system and depends on the molar fraction Xa- 
At fixed packing fraction, the larger the concentration of the big particles the larger is the viscosity. 

The agreement between numerical and analytical data is remarkable. At large packing fraction, a degree of departure 
from the analytical data indicates the major role played by the collisional integral that apparently is not fully captured 
by the numerical quadratures. This is an effect of the discretization and of the linear interpolation employed in order 
to compute some off-lattice points used in the calculation of the surface integrals. From the viewpoint of numerical 
stability, the range of packing fraction handled in simulation shows that the method is stable up to packing fractions 
of £3 < 0.4, rendering the approach usable in many experimental condensed matter contexts. In proximity of hard 
walls and under strong flow conditions, the numerical stability range reduces to packing fraction below £3 ~ 0.3. 



IV. FLOW IN A SLIT-LIKE PORE. TOY MODEL. 



We now consider a binary mixture of particles of equal masses and different sizes (a aa < &bb) confined in a narrow 
slit-like pore represented by two parallel smooth plates having infinite area and separated by a distance H + a aa (see 
fig. [2]). The walls are parallel to the yz plane and a flow along the z direction is induced by the presence of a uniform 
external field, F, also directed along z. Since H is the typical size of the system, u the average flow velocity and r\ 
the viscosity, the Reynolds number is Re = puH/rj. In the small Re regime, flows have negligible inertia forces. The 
viscous and pressure forces should be in approximate balance and a reduction, similar to the one occurring between 
the macroscopic Navier-Stokcs and the Stokes equation, takes place also at the present microscopic level of description. 

The first remark is that, besides being small in systems of microscopic size and in the limit of low flow velocities, 
the non linear term in the momentum equation vanishes. In fact, if one imposes the boundary conditions typical of a 
Poiseuillle flow, in a straight channel, it vanishes due to symmetry. 

We rewrite the l.h.s. of eq. (|T0| as: 



d t [p a (r , t)u« (r, t)} + Wi (p a (r , t)uf (r, t)u« (r, t) - p a (r , t)wf (r , t)w« (r , t)) 



vS (r , t) d tP a (r , t ) + V, (p a (r , i)< (r, t)) 



9 4 < + <(r,i)V i <(r,t) -V,[p a w lW 
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The first parenthesis in the r.h.s. is zero due to the continuity equation of species a, the second parenthesis is also 
zero due to the geometry of the problem and the stationarity of the flux, and the last term vanishes because w is 
directed along z and p a varies only along x. 

The flow induced by the presence of uniform fields parallel to the walls is characterized by velocities along the z 
direction u Q = (0,0,m"(x)) obeying the following set of stationary equations for each species: 



J d^(r,0[uf(0-ufM^ 



(48) 



In order to gain insight into the problem, we first set up and solve a toy model, obtained from the full model under 
some crude approximation. We replace the true non-uniform density profiles by two slabs of constant densities, 
nA,riB, respectively: 

n A {x) = n A e{x)6{H-x) (49) 
n B (x) = ub0{x — a)6(H — a — x). 

(50) 

As shown in ref. [36| the kinetic contribution to the viscosity due to the individual species is rj^J = ^rr n Ai while 
the collisional part is given by formula (|4"Tj) . Notice that the above prescription for the kinetic contribution to the 
viscosity is a consequence of our single time relaxation ansatz. Let us introduce the abbreviations: 

Vaa = (v a K a + V { aa) . Vbb = (vbb + Vbb) > VA = (vaa + Vba) > = (vbb + Vab)- ( 51 ) 

When the two species have different diameters, the smaller species can approach the wall up to a distance (7,4,4/2, 
whereas the larger B species can only reach the distance a AA /2 + a, where a = (jjbb — &aa)/2- In the following, we 
take H > 2a and measure distances, x, from the position of closest approach of particles A to the wall as illustrated 
in Fig. [2] 

We now determine the velocity profiles consistent with our toy model. In the central region a < x < (H — a) we 
obtain the following coupled differential equations for the partial velocities: 

VAA d^(x) + r] AB d*uf(x) = r[uf(x) - uf(x)} - F z A n A (52) 

TjBBdlufix) + r) AB d%u A (x) = -T[u A {x) - uf(x)] - Ff n B (53) 

with 

r = n A 7 AB = I ^2nK B Tii A B(T 2 A B9ABn A n B . (54) 

In the excluded layer for the big spheres, that is for x < a or x > (H — a), the equations read: 

rj AA d 2 x u A {x) = -F A n A , uf{x) = 0. (55) 

The solution with boundary conditions of zero velocity u A at x = and x = H 1 together with zero velocity uf at 
x = a and x — H — a, and continuous u A at x = a, reads: 

A . . F A n A (H~x)x , TT . ._„. 

u A {x) = — — , x < a, x>(H-a) (56) 

Vaa 2 

while for a < x < (H — a) is given by: 

4/ \ F A n A + F B n B AH — x)x {H — a)a, _ , , , w „,„^ . „,„„> 
= — . X ~b (" 1T~ - ~ ^)+^(cosh(A(.T - H/2)) - cosh(A(a - H/2))) 



and 



r)A + V 
F A n A (H - a)a 
Vaa 2 



= F * AnA ± % nB fiJL^h (JL^\+C B (cosh(X(x H/2)) cosh(A(a - H/2))) (58) 



(57) 



T]A + 1] 
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with 



r( ^ + 7 J (59) 



and 

t}b ( 1 F^n A f]B - Ffn B f]A F A n A (H ~ a)a^ 



A fj A +fj B \T va + Vb ' Vaa ' 2 J cosh(A(a - H/2)) 1 ' 



C B = -4^C A . (61) 

11B 

The main result of this illustrative model is the presence of different velocity profiles for the two species. The 
velocity profile of the bigger particles is generally lower than the corresponding profile of the smaller particles, because 
B effectively sees a narrower slit, so that the confinement effect resulting in a Poiseuille parabolic-like profile is more 
relevant than for the other species. The velocity difference between species remains finite in the center of the pore. 
Fig. [3] displays the velocity profiles of the two components for two different values of the composition and at low 
packing fraction. We will see below that even such a simplified model is able to capture some important features 
which are better studied by means of our full numerical method. 



V. LBM NUMERICAL RESULTS 



In the present section we solve numerically the equations for the phase space distributions and study the associated 
steady state properties. The most interesting applications of our method are those which regard truly inhomogeneous 
situations where one observes the interplay between the microscopic structure and the flow properties. The small 
amount of material in the confined region makes experimental probing of confined fluid extremely difficult and thus 
numerical methods are welcome. 



A. Profiles. 



We reconsider the Poiseuille problem of the previous section using the full power of our theory to obtain the density 
and velocity profiles of the two components when a uniform field parallel to the plates is applied. 

The available region for each hard sphere can be computed according to the density profiles that show where the 
centers of the spheres can range. Given that the sphere radii are multiples of the mesh spacing, the extrema of the 
available region can only fall on a mesh point. In this respect, the values of H are equal to 7, 13 and 29 mesh-points for 
the three channels, that is H/cfaa = I -75, 3.25 and 7.25 respectively. On the other hand, for the velocity profiles, we 
have added a mesh spacing on the left and right extrema to account for the non-slip boundary condition as imposed 
via the bounce-back scheme in the numerical method. When looking at the intervals from the velocity point of view, 
they look wider. 

Fig. 2] displays results for the density and velocity profiles for each species in a channel of width H = 1.75ctaa 
for a mixture of particles with diameters gaa = 4 and gbb = 8 in two different cases: packing fraction £3 = 0.084 
and composition Xa = 0.25 and packing fraction £3 = 0.073 and composition xa — 0.75. In the upper left panel 
the mixture is almost a gas and the density profiles for both compositions considered do not display any significant 
structure, whereas in the lower left panel the velocity profiles have nearly parabolic shapes, each species having its 
own curvature. One may conclude that in such a regime the flow behavior of the confined fluid is Poiseuille-like. As 
we increase the packing fraction to £3 = 0.42 (composition Xa = 0.25) and packing fraction £3 = 0.26 (composition 
Xa = 0.75), we observe the behavior reported in the right panels of Fig. 01 The large species B is unable to display 
significant peaks, due to the small plate separation that inhibits the formation of two layers of B particles. The small 
species, instead, starts developing some density enhancement in the proximity of the walls (upper right panel). The 
velocity profiles instead still bear a parabolic behavior with the smaller species always being faster than the bigger 
species. 

In Fig. [5] we consider the same state conditions but a wider slit, H = 3.25a aa, such that more than one layer of 
particles B can fill the gap. In the low packing region, at £3 = 0.073 and at £3 = 0.084 (upper left panel), both density 
profiles are rather flat and the associated velocity profiles do not display significant structure, the smaller particles 
being faster (lower left panel). As we consider larger packing fractions, £3 = 0.42 and composition Xa = 0.25, we 
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observe a stronger structure in the large particles than in the small particles. The other packing fraction, £3 = 0.26 
(composition Xa = 0.75), also produces some significant peak structure at the walls for the large particles, but the 
small particles have a rather flat density profile. The small rlS are somehow enslaved by the large component as 
far as the structure is concerned. In the velocity profiles (lower right panel) we observe that the B component for 
Xa = 0.25 has a higher velocity than the small counterpart near the center-line of the slit. This effect is determined 
by packing: small particles have more room near the walls while large particles move faster near the center-line. 

In Fig. [5] we consider a channel with H = 7.25cfaa- In the low packing region (left panels), there is not significant 
structure, the velocity profiles are Poiseuille-likc and display the same trend as in the corresponding panel of Figs. 
[4] and [5] More interesting is the situation at higher packing, shown in the upper right panel. When £3 = 0.42 
and Xa = 0.25 the large particles display four well defined peaks, while the small particles have a flatter structure. 
Correspondingly, the velocities of the two species (lower right panel) are anti-correlated with the density profiles, 
the largest velocities being attained where the density displays local minima. Also the velocity profiles at such large 
packing fraction show an inversion, that is the larger species has the larger velocity. The other case considered for 
£3 = 0.26, Xa = 0.75 displays a less pronounced structure and nearly no oscillations in the velocity profiles. 

To conclude, the density profiles are rather sensitive to the composition and show the characteristic enhancement 
near the walls as packing increases. On the contrary, we have not found a sensitive dependence of the density profiles 
on the applied load. The velocity profiles for moderate packing fractions have shapes reminiscent of the parabolic 
Poiseuille-like profiles, with the velocity of the larger species being smaller at low packing , in agreement with the 
prediction of the toy model of the previous section. 



B. Coarse grained observables. 



Besides considering the microscopic aspect represented by the various profiles, it is of interest to analyze some 
average properties, such as the volumetric flow of each component and the selectivity, quantities which are more 
easily measurable in experiments. We define the volumetric flow rate of each species as: 

Q a = Af dxn a {x)u a {x) (62) 
Jo 

where A is the area of the plates. In the upper panels of Figs. [7] and [8] we display how these quantities vary with 
the packing fraction, £3, for a fixed value of the mixture composition, Xa- In the lower panel we display how the 
selectivity, S, defined as 

Ss §W (63) 

varies with the packing at fixed composition. 

We first observe that the volumetric flows at low densities increase almost linearly with packing . In general, the 
volumetric flow of the small species for equal composition (Xa = 0.5) is larger for the small species, on account of the 
larger effective volume available. If packing increases further, the volumetric flow decreases. Such a feature has the 
following explanation: in the low density region, Q a increases almost linearly with packing, but as £3 becomes larger 
it determines a quadratic increase in viscosity and thus a decrease of the velocity of the fluid in the pore. As a result , 
there exists an optimal value of the packing for which the volumetric flow is maximum. The peak shifts towards 
larger values of the packing as the pore becomes larger. Interestingly, such a non-monotonic behavior is captured also 
by the toy model of section [TV] . The results suggest that particles could be separated by their size. Hydrodynamic 
chromatography, instead, exploits the Poiseuille velocity profile in channels Q- 



VI. CONCLUSIONS AND PERSPECTIVES 



To summarize, we have proposed a new theoretical method to study mixtures under confinement by using concepts 
derived from both DFT and LBM. The method has been applied to recover bulk properties such as the shear viscosity, 
and to study inhomogeneous situations, such as the flow in a slit-like channel of molecular thickness. Whereas MD 
simulations provides a complete microscopic picture for molecular level mechanisms at the price of a demanding 
computer time, the present method offers an efficient alternative in terms of CPU-time and high numerical versatility. 

We have extended the kinetic model described in refs. |3"il [i^ to a multicomponent system. The formalism enables 
us to compute not only the transport coefficients of fluid mixtures, but also to investigate inhomogeneous dynamical 
properties. We have verified that the numerical solution is in agreement with the values of viscosity predicted by our 
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theory. In addition, we have considered the dynamics of a mixture of particles of different size confined between two 
parallel plates. Under high confinement, different species display different velocity profiles. We also found non-trivial 
behavior when the plate separation becomes so narrow that the harshly repulsive forces among the particles and with 
the confining walls come into play. 

We wish to comment that our approach relies on the validity of the local equilibrium approximation and is cor- 
roborated by the observation that both local thermodynamic and hydrodynamic theories appear to be reliable up 
to surprisingly small wavelengths. At smaller scales, both conditions break down and a transition to a genuinely 
molecular regime takes place. Compared to previous approaches (40l - [4^ |. our method contains two main ideas which 
in our opinion render it a workable scheme: the simplified kinetic kernel (|15|) and the numerical scheme based on 
the LBM techniques. Very fast numerical solutions can be obtained along these lines and complex geometries can be 
handled with great simplicity and good analytical insight. 

In developing the numerical algorithm via a generalized version of the Lattice Boltzmann method, we did not take 
into account temperature variations across the system because we wanted to maintain the methodology as simple 
as possible. The non isothermal extension of the algorithm is doable, but this would require to include Hermite 
polynomials up to fourth order (that is, 127 mesh neighbors in the collisional part of the algorithm), in order to 
guarantee numerical stability and physical accuracy. In the future, we plan to handle this aspect. 

Before concluding we wish to add that several applications of our theory can be envisaged to the study of fluids in 
nanospaces. Among others we mention just a few: 

a) the spreading of a highly concentrated region in flow conditions and of the competition between the advection 
and molecular diffusion [531 ]: 

b) the transport of different molecular species in a channel of varying width and the influence of the size and/or 
mass on the dynamical properties Q; 

c) the entropic, Asakura-Osawa interactions with walls and the interplay between advection and diffusion. When 
particles of different size are present one expects that the larger particles experience a larger attraction towards the 
walls, the so called effective entropic force 0, [55[. At low Reynolds number, the scenario should not be greatly 
modified, but the issue is worth to be considered with some care; 

d) the dynamics of menisci formed in narrow slits between two coexisting liquid phases (60j ; 

e) driven flows past chemically, physically heterogeneous substrates, a mechanism which could enhance mixing [57j : 

f) since the equilibrium phase diagram of substances under confinement is shifted with respect to bulk conditions, 
a similar dynamical effect is expected and could be used to obtain specific properties [58j | ; 

g) cylindrical pores are known to lead to strong, non-trivial size selectivity, which is highly sensitive to the pore 
width [59| . It would be very interesting to exploit our model in such a geometry. 



Appendix A: The connection with the chemical potential 

In this appendix we give a proof of formula (j3"5")) . To this purpose let us consider the r.h.s. of the momentum 
equation (fTU|) for species a. At equilibrium the l.h.s. of eq. (fTUf vanishes together with the drag and the viscous 
contributions, so that we obtain the condition: 

k B TVn a (r, t) - n a (r, t)F a (r, t) = n a (r, t)F a ' mf (r, t) (Al) 

Under the adiabatic hypothesis, we assume the non-equilibrium positional pair correlations to be given by the 
equilibrium pair correlations [2(1 [23] evaluated when the instantaneous density profiles assume the values n A (r, t) and 
n B (r,t). The force F a ' m '(r) can be expressed as 

F«.™/(r) = - / dr'^. 9 ^(r,r')^(r')Vt/^(|r-r'|) (A2) 
J P 

where we have used the property of the hard sphere potential: 

j dr'nP(r')g a p(r,r>)V r U^(\r-r'\) = -k B T J dr'n (r')g a0 (r, v')6(\v - r'| - ^O^^r (A3) 

in order to recast the surface integral over the sphere of radius a a p eq. (|35[) into a volume integral. Substituting the 
result (|A1[) into eq. (|A2[) we obtain the Born- Green- Yvon relation (25j. We now compare the r.h.s. of eq. (|A1[) with 
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the exact equilibrium relation 




where Tint is the non-ideal contribution to the free energy functional, — fcsT 1 c Q( g(r, r') = sn^^^TP) * s ^ e Ornstein- 
Zcrnikc direct correlation function and — V/^ 4 is the force arising from the interactions with the other fluid particles 
on a particle located at r [43j]. By equating the r.h.s. terms of eqs. (|A2[) and (|A4[ we conclude that eq. (f3"5)) holds. 
Finally, notice that eq. (|A4[) can be rewritten as 

fc s TVlnn«(r) + V/x?„ t (r) = F a (r) (A5) 

that is, in thermal and mechanical equilibrium the net force on a fluid particle at position r must vanish. This occurs 
when the entropic force, represented by the first term, being the resultant of the attractive and repulsive forces exerted 
on a molecule at r from all the other molecules (the second term on the l.h.s.) and the external force (on the r.h.s.) 
compensate exactly. 



Appendix B: The discretization procedure 



A crucial goal of the present work is to achieve the solution of the coupled integro-differential equations for non- 
ideal mixtures by means of a suitable numerical scheme. In this work, we consider the Lattice Boltzmann method as 
the reference framework that enables us to solve the kinetic equation in bulk conditions and under confinement, in 
particular when the confining geometry is highly non-trivial. The main asset of the LBM is to work directly with the 
distribution functions by accomplishing a spectral decomposition of the distribution in velocity space via a Hcrmitc 
representation. In this way, the three-dimensional velocity space is reduced to a handful of sampling points and the 
three-dimensional real space is discretized over a Cartesian mesh. The resulting distribution function is transformed 
into a reduced set of populations. The phase space evolution is thus rewritten as an efficient updating algorithm 
where the streaming operator has the form of a forward Euler updating step. The LBM approach has been previously 
employed to solve the collisional dynamics of the single-component hard-sphere fluid and provided robust solutions 
up to packing fraction < 0.35 [3tjj . 

The LBM is a very general framework to handle the evolution of generic kinetic equations and does not depend in 
any specific way on the form of the collisional kernel [lj| . We exploit this generality to solve the collisional dynamics 
encoded by eq. (|15[) . where non-local structural forces appear as surface integrals, and retain the simple LBM picture 
of updating the populations. In addition, the BGK kernel of cq. (|15p also depends on the Hcrmitc representation 
via the equilibrium distribution. A second-order expansion of the distribution function, corresponding to a Hcrmitc 
representation to second order, is usually employed in one-component LBM implementations [38l I6ll463| . 

The distribution function of each species is discretized in velocity space and the continuous phase space distribution 
is replaced by a Q-dimensional array, / Q (r, v, i) — > f p (r,t) where p = 1,...,Q labels discrete velocities c p . The Q 
discrete velocities connect neighboring mesh points on a lattice and mirror the hop of particles between mesh points, 
generally augmented by a null vector Co accounting for particles at rest. The form of the lattice velocities c p depends 
on the order of accuracy of the method and reflects into the required Hermite order [frij . 

In essence, the LBM dynamics is achieved by a rearrangement of populations over spatial shifts via an explicit 
update arising from the following approximation to the streaming operator, 

w ,, tvW ,, B itia*^L4M (B1) 

where At is the time-step. We choose the mesh classified as D3Q19, that is, containing 19 discrete velocities connecting 
first and second neighbors of a cubic lattice and overall achieving second-order Hermite accuracy. We choose equal 
masses for the hard-sphere species so that the thermal velocity is the same for both species and given by Vt = I/a/3- 
An extension of the present method to different masses can be easily obtained by using the method described in ref. 
[221 ] . The updating step for the populations is compactly written as 



(B2) 
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where /™(r, t) is the post-collisional population of species a, containing the effect of the BGK relaxation and the 
hard-core collisions, written as 



The BGK contribution to the post-collisional population has the form 

f^ a (r,t) = (1 - W At)/ p a (r,t) +^A^(r,i) 

where 

Cpiuf (r, t) (cpiCpj ~ v^Sij) (v^Sij + 2uf(r, t)uj(r, t) - u t (r 7 t)uj(r, t)) 



^lJr,t)=w p n a (r,t) 



1 
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and 



^(r,t)=w p n a (r,t) 



CpiUi(r,t) (cp^ - i>jAj)uj(r,f)uj(r,f) 



2v\ 



(B3) 
(B4) 

(B5) 
(B6) 



The w p are the Q weights normalized to unity that arise from the Hermite expansion [19| . Formula (|B6[) is a low-Mach 
(0[A/a 3 ]) expansion of the local Maxwcllian corresponding to the barycentric velocity for a mixture of equal masses 



u(r,t) 



£ Q n a u a (r,t) 



with 



n a (r,t) = J2J pM) 
p 

* a (r,t)u a (r,t) = ^c p / p «(r,t) 



(B7) 



(B8) 



The hard-sphere post-collisional contribution is 

, pi c[ a) {v,t) , (c 



# a,a (r,t) 



4%) Ui (r,t)C a) (r,t) 



(B9) 



The collisional integral is approximated via the following real-space quadrature 

C (Q) (r, t) = J2 C al3 (r, t) = -4tt ^ a% ^ W,jk 9ffQ/ j(r + |a Q ^k g , t)n„(r, t)np(r + <J a pk q , t) 



4 _ 2^/^k 9 • [^(r + a^k,,*) - u«(r,t)] 



(BIO) 



where k 9 are the nodes of the quadrature over a spherical surface [65| and W q the associated quadrature weights. 
We choose a 18-points quadrature that is exact for quadratic integrands. The quadrature nodes on the unit sphere 
are the six on-axis points (±1,0,0), (0,±1,0) and (0,0, ±1) and twelve points along the diagonals (±-\/2, ±v2, 0), 
(±\/2, 0, ±\/2) and (0, ±y/2, ±y/2). Correspondingly, the quadrature weights are W q = 1/15 for the on-axis points 
and W q — 1/30 for the diagonal points. It should be borne in mind that the quadrature points do not necessarily 
fall on the mesh. For instance, if the the HS radius is an integer and one considers intraspecies collisions, there are 
only six quadrature nodes that fall on the mesh while twelve nodes along the diagonals do not. In order to sample 
the fluid density and velocity on off- mesh points, we employ a trilinear interpolation scheme (66j . It is clear that 
the accuracy of the collisional integrals improve for larger HS diameters, since the eight mesh points employed in the 
trilinear interpolation fall closer to the spherical surface. The same type of consideration applies for the interspecies 
integral and, in all circumstances, for the midpoint rule employed to estimate the radial distribution function via the 
Carnahan-Starling expression eq. (|44[) . In our numerical applications, we applied a trilinear interpolation for any 
off-mesh quadrature nodes and midpoints for the radial distribution function, and chose HS diameters to be even 
numbers in order to minimize the number of off-mesh interpolations. 
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In solving eq. (fT5|) with the LBM approach, we need to control the numerical error introduced by the discretization 
procedure over velocity, space and time variables. It is well-known that for one-component, uncorrelated fluids LBM 
bears three types of error as compared to the analytical solution of the incompressible Navier-Stokes equation. The 
first error depends on the mesh spacing Ax and decays quadratically with Ax. In presence of hard boundaries, and 
in particular with curved fluid- wall interfaces that are approximated as staircase geometries, the accuracy degrades 
to linear dependence with the mesh spacing. The second type of error depends quadratically on the time-step At. 
Finally, the so-called compressibility error arises from the fact that LBM docs not exactly enforce a non-zero divergence 
of the velocity field, with a resulting error that depends quadratically on the Mach number. As a side effect, the 
time-step error also degrades linearly with At. These three types of error can be made arbitrarily small by changing 
Ax, At and the mass unit while keeping fixed the physical value of the mesh spacing, kinematic viscosity and mass 
density. In order to minimize the composite error, the mesh spacing is typically the parameter that is reduced, by 
negotiating with the numerical effort needed to resolve local morphological details and flow pattern. A concomitant 
effect of the space-time discretization is the effective reduction of the kinematic viscosity arising from a negative 
viscosity of numerical origin, — At ^ that depends on the specific form of the lattice employed, such as D3Q19 in our 
case. The numerical contribution to viscosity is evaluated via a Chapman-Enskog analysis with the outcome that the 
effective viscosity can be controlled to arbitrary accuracy [l^. For the case of correlated dynamics, the viscosity due 
to the finite-time propagation enters with the same numerical value of the uncorrelated case, as shown in our previous 
publication j36j . 

The extension of the LBM to non-ideal fluids, such as the binary mixture described in this paper, carries the same 
types of error as for the uncorrelated case. The BGK contribution to the mixture dynamics, encoded by the kernel 
cq. (|15|) . generates the same type of numerical error as for the single-component case. The issue of resolving the 
atomic-scale structural correlations critically depends on the rapidly varying oscillations of the density profiles, as 
much as on the current and higher moments profiles. The profiles need to be resolved at spatial scale smaller than 
the molecular radius Aa; < min Q a aa /2. When confronting with experiments performed on nanoscopic systems, the 
computational load needed to solve the kinetic equation requires a number of mesh points that scales as ~ <7^ n . 
The same argument applies when resolving the temporal evolution, with a time-step that must be smaller than the 
typical collision time, At < min a <7 aa /vT- Regarding the additional compressibility arising in LBM, the error plays a 
minor role as compared to the highly compressible nature of fluids at molecular scale. The low-Mach numbers often 
employed in nanofluidics, condensed-matter conditions further alleviates this problem. 

In addition to the standard error terms, the computation of the collisional integrals eq. (|B10[) introduces a further 
source of error. It is a simple exercise to show that for a uniform system the diagonal points in the quadratures 
produce the collisional contribution to viscosity. Therefore, the numerical error on the collisional viscosity are mostly 
due to the off-lattice location of the quadrature points and the employed trilinear interpolation scheme. The level of 
accuracy of viscosity can be controlled in a systematic way and the theoretical curve can be recovered to an excellent 
level. Our numerical results indicate that, once the molecular spatial scale is resolved to account for the microscopic 
flow patterns, the error due to the quadratures becomes negligible. 

Lattice nodes adjacent to the walls obey collision rules different from those characterizing bulk nodes. We have 
adopted the so called bounce-back collision rule which states that the velocity of a particle incident on a wall is 
reversed after the collision. 
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FIG. 1: Shear viscosity (normalized to the ideal gas value) versus packing fraction for a bulk mixture at different values of 
the composition Xa = n A (jiA + tlb) of hard spheres of diameter ratio <jbb/&aa = 2 . The viscosity decreases with the 
concentration Xa- The continuous curves represent the values obtained with the theoretical formula, whereas the symbols refer 
to the numerical values obtained from our LBM computer simulations. 
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FIG. 2: Sketchy view of the system. The walls, parallel to the yz plane are located at positions x 
the flow occurs along the vertical z direction. The direction y is normal to the figure. 



= -ctaa/2 and x = H+aAA/2, 
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FIG. 3: Velocity profiles of the two species for a channel of width H — 6<taa and load F = Q.QQlksT /gaa, according to the 
toy model. The ratio of the diameters of the spheres obbIoaa = 2. Velocity of species A (open black circles) and of species B 
(red squares) at bulk composition Xa = 0.25 and packing fraction £3 = 0.084. Velocity profile of species A (green diamonds) 
and of species B (blue triangles ) at bulk composition Xa = 0.75 and packing fraction £3 = 0.073. 
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FIG. 4: Numerical results for a channel of width H/cjaa = 1-75 and load F — W~ 6 kBT/aAA- The spheres ratio of diameters 
is obb/o-aa = 2. Upper left panel: Density profile of species A (open black circles) and of species B (red squares) at bulk 
composition Xa = 0.25 and packing fraction £3 = 0.084. Density profile of species A (green diamonds) and of species B (blue 
triangles ) at bulk composition Xa = 0.75 and packing fraction £3 = 0.073. Lower left panel: Velocity of species A (open black 
circles) and of species B (red squares) at bulk composition Xa = 0.25 and packing fraction £3 = 0.084. Velocity profile of 
species A (green diamonds) and of species B (blue triangles ) at bulk composition Xa = 0.75 and packing fraction £3 = 0.073. 
Upper right panel: Density profile of species A (open black circles) and and of species B (red squares) at bulk composition 
Xa = 0.25 and packing fraction £3 = 0.42. Density profile of species A (green diamonds) and of species B (blue triangles ) at 
bulk composition Xa — 0.75 and packing fraction £3 = 0.26. Lower right panel: Velocity of species A (open black circles) , of 
species B (red squares) at bulk composition Xa = 0.25 and packing fraction £3 = 0.42, velocity of species A (green diamonds) 
and of species B (blue triangles ) at bulk composition Xa = 0.75 and packing fraction £3 = 0.26. 




FIG. 5: Numerical results for a channel of width H/ctaa = 3.25 and load F = W~ 6 kBT/tJAA- The ratio of the diameters of 
the spheres is obbI^aa = 2 . The symbols used in the four panels correspond to the same values of packing and concentration 
as those employed in figure [4] 
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FIG. 7: Volumetric flow rates for species A (upper left panel) and B (upper right panel) and selectivity function S (lower panel) 
as a function of the packing fraction for a channel of width H — 2a aa and Xa = 0.5. For the sake of comparison we plotted 
the lines obtained assuming a simple Poiseuille type of behavior of the system as explained in the text. 
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FIG. 8: Volumetric flow rates for species A (upper left panel) and B (upper right panel) and selectivity function S (lower panel) 
as a function of the packing fraction for a channel of width H — 2a aa and Xa = 0.5. For the sake of comparison we plotted 
the lines obtained assuming a simple Poiseuille type of behavior of the system as explained in the text. 



